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ABSTRACT 



The nonlinear clustering of dark matter particles in an expanding universe is usually 
studied by N-body simulations. One can gain some insight into this complex problem 
if simple relations between physical quantities in the linear and nonlinear regimes can be 
extracted from the results of N-body simulations. Hamilton et al. (1991) and Nityananda 
and Padmanabhan (1994) have made an attempt in this direction by relating the mean 
relative pair velocities to the mean correlation function in a useful manner. We investigate 
this relation and other closely related issues in detail for the case of six different power 
spectra: power laws with spectral indexes n = —2, —1, cold dark matter (CDM), and hot 
dark matter models with density parameter = 1; CDM including a cosmological constant 
(A) with VLcdm = 0.4, A = 0.6; and n = -1 model with O = 0.1. We find that: (i) 
Power law spectra lead to self-similar evolution in an O = 1 universe, (ii) Stable clustering 
does not hold in an O = 1 universe to the extent our simulations can ascertain, (iii) Stable 
clustering is a better approximation in the case of O < 1 universe in which structure 
formation freezes out at some low redshift. (iv) The relation between dimensionless 
pair velocity and the mean correlation function, £, is only approximately independent 
of the shape of the power spectrum. At the nonlinear end, the asymptotic value of the 
dimensionless pair velocity decreases with increasing small scale power, because the stable 
clustering assumption is not universally true, (v) The relation between the evolved £ and 
the linear regime £ is also not universal but shows a weak spectrum dependence. We 
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present simple theoretical arguments for these conclusions. 

Cosmology: large-scale structure of Universe - cosmology: theory - galaxies: clustering - 
numerical methods 
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1. INTRODUCTION 

Large-scale structures (like galaxies, clusters etc.) are believed to have evolved primarily 
through the action of gravitational force out of small initial density perturbations. When 
these density perturbations are small, it is possible to study their evolution using linear 
theory. But once the density contrast becomes comparable to unity, linear perturbation 
theory breaks down and one must invoke either analytic approximations or direct N-body 
simulations to study the growth of perturbations. 

Several people have performed extensive numerical simulations, with varying dynamic 
range and resolution, to study this problem. While these simulations are of value in making 
concrete predictions for specific models, they do not seem to have provided deep insight into 
physics of non-linear gravitational dynamics. To obtain such an insight into this complex 
problem, it is important to ask whether simple patterns in gravitational clustering can 
be extracted from the plethora of simulations. If such patterns are found, they will work 
as pointers for possible analytic approximations and could even eliminate the need for 
extensive simulations in some specific cases. We shall study, in this paper, the possible 
existence of such simple patterns in the nonlinear gravitational clustering. 

It would be naive to hope for a simple picture to emerge in the full theory of structure 
formation, involving both gravitational and hydrodynamic processes. Radiative coupling 
of baryons and photons is a multiscale process with several dimensionless ratios, and it 
is unlikely that such processes can be described by simple universal rules. On scales less 
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than a few Mpc such physical effects become increasingly important. The situation is 
quite different as regards the dark matter in the universe, often assumed to be composed 
of collisionless, weakly interacting, massive particles. By limiting our consideration to the 
gravitational interaction alone, one may hope to uncover some basic relationships. 

Section 2 discusses some possible scaling relations and motivates the numerical 
simulations we have done. Section 3 describes the simulations and last section presents 
the results and discussions. 



For the case of interest, we can describe the dark matter as a system of particles interacting 
via Newtonian gravity, at scales which are small compared to the Hubble radius. The 
motion of particles in this case can be described by the equations 



where Xj (£) is the comoving position of the i-th particle, a (t) is the expansion factor, 
Pb (t) is the background density in a matter dominated universe and is the gravitational 
potential due to the perturbed density dp = p — pb = Pbd- All spatial derivatives are with 
respect to the comoving coordinates. In the fluid limit, the same system is described by 
the equations 
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Xi H x, = ^V0; V 2 (p = 4nGp b a 2 8 : 



a a 



(i) 



6 + V • [(1 + S) u] = 0; u + (u • V) u + — u = -^V0, 

a or 



(2) 
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where u = (dx/dt). For a universe with density parameter fi = 1, we can introduce the 
time coordinate r = In (t/T) where T is an arbitrary positive constant and a rescaled 
potential ip = (4/9) Hq 2 (a<j)) with H the Hubble constant. Then equation (1) reduces to 

; /-X 1 + lrfx 1 = _ v ^_ V 2 tf; = ^ (3) 



dr 2 3 dr 



while (2) becomes 



^ + V-[(l + 5)p]=0; |B + (p.V)p+ip = -Vik (4) 

where p is the velocity field corresponding to r coordinate. All reference to the background 
spacetime has completely disappeared in these equations. Hence these equations have no 
intrinsic scale. If we now further assume that the initial power spectrum of perturbations 
is scale invariant with P(k) = Ak n , then we expect the dynamical evolution of such a 
system to possess some universal characteristics. 

To avoid misunderstanding, we stress several points. First, any particular realization 
of an initial power spectrum will have scales appearing in it randomly. When averaged over 
an ensemble of realizations with the same P (k), these scales will disappear. Second, the use 
of a grid to generate initial conditions will introduce a scale of order the mean interparticle 
distance. Such scales are erased by evolution, except in the most underdense void regions. 
Third, softening in the gravitational force introduces a minimum resolution scale which 
cannot be avoided. We shall, however, concentrate on scales which are larger than 
the softening scale to obtain reliable conclusions. Fourth, requirements of convergence, 
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(n > —3), and the generation of a k 4 tail due to discreteness effects restrict the range of 
n to — 3 < n < 4. If equations (2) are studied perturbatively, spurious divergences due 
to long wavelength modes can arise for — 3 < n < — 1. The long wavelength contribution 
to bulk velocity, for example, has a divergent contribution for n = — 2 when handled 
perturbatively. The nonperturbative argument given in the paragraph above clearly shows 
that these divergences are spurious. Thus endpoint effects should not be important for 
— 3 < n < 4 and thus no physical scales enter the problem for Q = 1. 

One immediate consequence of the above analysis is that a simple scaling relation 
must exist in all simulations involving dark matter with power law spectrum, in an O = 1 
universe. Since there is no a priori scale in the problem, the natural length at any time is 
^nl(«), the comoving scale which is going nonlinear at the epoch a . Thus all dimensionless 
physical parameters must be expressible as a universal function of (cc/^nl)- For example, 
we expect the correlation function £,(x,a) to have the form a) = f n [x/x^(a)}. This 
implies that the evolution is self-similar, and £ is completely specified by a single function 
f n which depends on the index n. 

Further progress can be achieved if we assume that clustering at the extreme nonlinear 
end is "stable" (Peebles 1980). In such a case, £(a, x) must have the form ^(a, x) = a s F(ax) 
in the extreme nonlinear limit, because the number density of objects in collapsed regions 
is assumed constant, while the average density of the universe decreases as a -3 . This 
assumption is equivalent to the postulate that highly overdense, virialized objects maintain 
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their identities, i.e., merging and fragmentation either do not occur or occur at strictly 
balancing rates. The assumption of self-similarity requires that the correlation function at 
the linear end evolves as £l oc a 2 x~( n+3 \ Combining this result with the stable clustering 
postulate, it is easy to show that 

f (a, x) OC a 3 (ax)- 7 OC a 6/(n+5) x -3(n+3)/(n+5) ? ^ 

in the extreme nonlinear end. In other words, the index of the nonlinear correlation 
function 7 should be related to the index of the linear power spectrum n by 

3(n + 3) (6 ) 
7 (n + 5)' W 

if the stable clustering ansatz were true. Therefore, given the correlation function in the 
linear theory, it is possible to determine the correlation function in the nonlinear theory 
up to a constant factor provided we assume the validity of stable clustering. In section 4 
we examine the quantitative accuracy of the above two assumptions. 

Equation (6), however, is not very useful in practice. To begin with, it works only in 
the extreme nonlinear regime and one is often interested in the intermediate scales. Second, 
it does not provide any clue to the amplitude of the nonlinear correlation function in terms 
of the linear one. Third, it is based on the assumption that the power spectrum is strictly 
scale free. And finally, it rests on the ad hoc assumption of stable clustering. 

Attempts have been made to generalize this idea by using an effective index n e ff as a 
function of scale, defined as, 

dlnP 
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Nityananda & Padmanabhan (1994, hereafter NP) found that a simple implementation 
of scaling using n e ff does not work satisfactorily in that a smooth transition between 
linear and non-linear regimes is not obtained. It is important, therefore, to look for more 
sophisticated - yet manageable - rules which will provide information about the nonlinear 
theory in terms of the linear theory. 

One such attempt was made in NP along the following lines: They assumed that the 
mean relative pair velocity of particles v p (a,x) can be expressed in the form 

v p (a,x) = —axh[£(a,x)], (8) 

where £ is the average correlation function within a sphere of radius x and h is a universal 
function of its argument. Given this ansatz, NP show that it is possible to relate the mean 
correlation function £ to the mean linear correlation function £x, by the relation 

-£(«,*) dv 



(2 fl 
£ L (a,Z) = exp I - / 



(9) 



h(u)(l + u) J ' 

with / = x[l + ^(a,^)] 1 / 3 . Earlier work by Hamilton et al. (1991) suggested that such a 
relationship seems to be borne out by the numerical simulations. In other words, numerical 
simulations seem to suggest that the function h(£) is universal; it is independent of the 
epoch and the power spectrum used in the simulations. This function seems to behave as 
follows: h (£) (2/3) £ for £ <C 1, reaches a maximum of about 2 at £ ps 20 and decreases 
to about unity at large £. For a wide class of power law spectra and spectra similar to 
that of cold dark matter (CDM), the relationship in (5) can also be expressed as 
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f£ L (a,Z) (for £ L < 1.2, £< 1.2) 

£(a, x) = < 0.7[&(a, 0] 3 ( for 1-2 < £l < 6.5, L2 < £ < 195) . (10) 
1 11.7[£ L (a, /)] 3 / 2 (for 6.5 < £ L , 195 < £) 

This fitting function is due to Bagla and Padmanabhan (1993) and captures the essence 
of a more complicated fit given in Hamilton et al (1991). 

Unlike the self-similarity argument mentioned earlier, this relation cannot be 
independent of the power spectrum in the strict sense. To see this, consider a power 
spectrum with a sharp maximum at x ~ L and very little small scale power. [This is 
similar to the HDM spectrum]. The first scales which go nonlinear will correspond to 
clusters with L. The evolution of the first nonlinear structures in such a model 

could be well approximated by collapse of a spherical top hat (STH), since there is very 
little small scale power. For an STH model, one can easily show that h is a monotonically 
increasing function of the density contrast 5 ~ £ and h ~ £ 1,/2 as £ — > oo (see the discussion 
in section 4 below) . This suggests that h will increase with £ for a good range of £ in the 
HDM-like models until significant amount of small scale power develops. One can reach the 
same conclusion from (10) as well: If there is very little small scale power in linear theory, 
then £ will have to rise far more steeply with £^ (in the nonlinear end) than suggested by 
(10). 

The study of such an extreme example suggests that there might exist a weak spectrum 
dependence in the h (£) relation even in more moderate cases. Broadly speaking, we expect 
h to be lower at a given £ (in the nonlinear regime with £ > 10 or so ) as we add more small 
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scale power. If stable clustering is invoked, then it is possible to motivate this conclusion 
along the following lines: 

Consider a spherical region of initial radius and overdensity vo where a = 
aor~^ n+3 ^ 2 is the variance of the gaussian density fluid. In a spherical model, this 
region will expand to a maximum radius of about (rj/W) and virialise to a final radius 
r = X(ri/ua) = (X/uao)r i 2 where A w 0.5. We shall assume that the coorelation function 
at the nonlinear end $,nl(^) is contributed by such virialised objects and can be computed 
as 

l+£(r)^(r) = <(^) 3 >. (11) 

This assumption is equivalent to the stable clustering hypothesis. From r = 
(A/z/aoK (n+5)/2 , it follows that 

f( r ) = ^-£J r -T^W ^6/(n+5)\ 

Assuming that v is a gaussian variable, we get 

Equations (12) and (13) show that £(r) oc £^/ 2 [/] with / pa r^ 1 / 3 ; however, the 
proportionality constant has a weak n-dependence. 
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(12) 



The averaging in (11) can be made more sophisticated by using a weightage 
proportional to r™. In that case, we still obtain the same r dependence but the gamma 
function becomes 



This reduces to the above result when m = 0; recently, Mo et al. (1995) have suggested a 
model based on m = 3, which - of course - leads to similar conclusions. 

While the above arguments suggest that h[£\ will have a weak spectrum dependence, 
they depend crucially on the assumption of stable clustering. The £ (r) can be modelled 
by virialised clusters only if stable clustering is valid. This statement is independent of 
the actual averaging process used in (12) . Because of this reason, arguments like the 
above ones - and the n dependence obtained from them - are unreliable. If the system 
does not obey stable clustering criterion, h[£] may still show an n- dependence but it cannot 
be obtained from the above argument. On the other hand, a theoretical explanation for 
the approximate universality (and the deviations from it) will go a long way in providing 
an analytic description of the nonlinear evolution of gravitating particles. It is, therefore, 
important to test this ansatz as precisely as possible in as diverse cases as possible. In 
particular, one must ask: (i) Do the relations h = h (£) and equations (10) have a weak 
spectrum dependence for pure power law spectra in O = 1 universe? (ii) How are the 
relations modified if the spectrum is not a pure power law? and (iii) What happens when 




(14) 
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the background universe is not described by a matter dominated, = 1 model? We 
present the results of such a numerical investigation in this paper. 

3. NUMERICAL MODELING 

Two kinds of N-body simulations are used to simulate the evolution of the particles. 
The first one uses a standard Particles- Mesh code (PM, cf. Hockney & Eastwood 1981; 
Efstathiou et al. 1985) with the staggered mesh scheme (see, e.g., Park 1990; Cen 1992). 
We use 240 3 = 10 7 1 particles on a 720 3 mesh with a simulation box size of 128/i _1 Mpc 
giving a spatial resolution of 0.18/i _1 Mpc and mass resolution of 4.2 x lO lo O/i _1 M . 
Gravitational force is calculated using an FFT technique coupled with the cloud-in-cell 
scheme to calculate the density as well as to assign the gravitational force to the positions 
of particles. 

The second kind of simulations uses the P3MG3A code (Brieu, Summers, & Ostriker 
1995). This code is based on the Particle-Particle-Particle-Mesh (P 3 M) scheme (Efstathiou 
and Eastwood 1981; Efstathiou et al. 1985) utilizing a special gravity processor called 
GRAPE 3A (Okumura et al. 1993) to compute the PP part. Each simulation uses 
128 3 = 10 6 - 3 particles on a 256 3 mesh with a simulation box size of 128/i _1 Mpc. The 
Plummer softening length is set to be 1/20 of the mesh cell size, giving a nominal spatial 
resolution of 0.025/i _1 Mpc and a mass resolution of 2.8 x 1O 11 O/i _1 M . 

Six different models are considered and are listed in Table 1. In the table, L box refers to 
the one dimensional size of the computational cube, l res lists the nominal resolution scale, 
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and as gives the linear value of the rms fluctuations in the mass on a sphere of radius 
8 h~ x Mpc at redshift zero. We adopt the transfer functions of Bardeen et al. (1986) for 
the standard CDM and HDM models. The power spectrum of the CDM+A model is 
from Cen, Gnedin, & Ostriker (1993). The initial density fields for each simulation are 
generated assuming Gaussian fluctuations and the initial velocity field is obtained using 
the Zel'dovich approximation. 

4. RESULTS AND DISCUSSION 

The above simulations were utilized to both examine the assumptions and explore 
the consequences of the arguments in §2. The important points to address are the self- 
similarity of scale free simulations, the validity of stable clustering, the universal nature 
of the relation between the dimensionless pair velocity and the mean correlation function, 
and the dependencies on spectral index. We discuss each of these features below: 
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4-1 Self-similar evolution 



We begin by verifying the self-similar nature of the evolution for power law spectra by 
testing whether £(a,z) can be expressed as £,[x/x n i (a)]. This question was investigated 
earlier by Efstathiou et al. (1988) who obtained results which were somewhat different from 
ours. We have now repeated this analysis with better resolution and range for n = — 2, — 1 
and we find that self- similarity is borne out in our simulations as expected. For the case 
n = — 1 the fit is excellent. For n = — 2 it is less good but, we believe, adequate. Figures 
(la,b) show the two-point correlation function as function of x/x n i (where x n \ is the scale 
corresponding to £ = 1 at three different redshifts z = 0.0, 0.5, 1.0 for n = —2, —1 models. 
Note that the curves fall on top of each other. Also shown (as two dashed lines in each 
figure) are asymptotic slopes at the extreme nonlinear and linear ends. This result is 
expected based on the theoretical arguments given in section 2, but serves as a good check 
on the accuracy of simulations. In addition, it clarifies certain doubts raised in Efstathiou 
et al. (1988) regarding self-similarity for the case n = —2. On all the figures we only 
plot the points which correspond to scales (separations) larger than 2.5 cells in the PM 
simulations and twice the softening lengths in the P 3 M simulations, which represent the 
true resolutions of two types of simulations. 

It may be noted that in the case of n = — 1, the linear slope extends to the quasi-linear 
regime as well, up to about £ ~ 10. This fact can be understood from equation (10). In 
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the quasi-linear regime, this equation predicts the evolution 

6 3(n+3) 

£oza~^+Vx l^W, (15) 

so that £ oc a 2 x~ 2 for n = — 1 in the quasi-linear regime as well. In other words, the linear 
and quasi-linear evolution has the same a and x dependence for n = — 1. (See Bagla and 
Padmanabhan, 1995 for a more detailed discussion of this point.) 

For comparison, We also show in Figures (lc,d,e,f) the corresponding results for the 
other four models (listed in Table 1): SCDM (lc), HDM (Id), LCDM (le), and open 
Pk = k~ x (If) models, respectively. All these models have a characteristic scale and the 
low-fi models have that feature even at the highest redshift probed. This should lead 
to deviations from self-similarity in figures 1 c - e. Such a deviation is apparent in the 
case of HDM (which has the strongest scale dependence), somewhat less apparent but 
noticeable in the case of CDM + A and CDM (in which the scale dependence is weaker) 
and ambiguous in the case of open n = — 1. This is essentially due to the fact mentioned 
above: n = — 1 is special in the sense that linear evolution holds even for the quasi-linear 
regime. Hence the deviations are somewhat offset in this case. 

In a truly self-similar evolution, we also expect velocities of particles to be a function 
of [x/x n i (a)]. To test this we studied the behavior of the ratio i>3d(a, x)/[ax non u near (a)] 
of the 3-dimensional velocity dispersion to our unit of velocity. Figures (2a,b,c,d,e,f) show 
the results for the different cases. It is seen that self-similarity (and the deviations from 
it) has the same pattern for velocity dispersion as for correlation functions: Figures (2a,b) 
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show time invariant behavior but (2c,d,e,f) do not. Note that Figure 2f shows the non- 
self-similar behavior of the open n = — 1 model that is hidden in Figure If. Note also the 
somewhat surprising result that in all cases presented, the velocity difference between two 
particles is a rather weak function of their separation within 10 ±15 of the nonlinear scale, 
especially for the self-similar cases. 

4-2 Relative pair velocity, 'universality' and stable clustering 

According to the arguments presented in section 2, we expect curve h(£) to vary slightly 
but systematically as more small scale power is added. We indeed find such a behavior. 

Figures (3a,b) show the behavior of h for n = —2 and n = — 1 at three different 
redshifts z = 1, 0.5, 0. In the nonlinear end, the curves for n = — 1 is lower than those for 
n = —2. That is, as we add more small scale power, the relative pair velocity decreases at 
nonlinear scales. This behavior is further borne out in the case of CDM (see Figure 3c) 
and HDM (see Figure 3d). In the redshift range considered, CDM model has an effective 
index which is between n=-l and n=-2 and HDM has less small scale power than n = —2. 
The curve for CDM falls in between those for n = — 1 and n = — 2 and that of HDM is 
above the one for n = —2. 

The behavior of /i(£) for the HDM model allows a simple test of the reasoning presented 
in section 2. To do this properly let us consider the evolution of a spherical top hat model 
in a background universe with some O along the lines of Gunn and Gott (1972). When 
O 7^ 1, we have a (t) oc (l — O -1 ) where Q now refers to the instantaneous value at time t. 
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The quantity h can be estimated in the STH model by the relation (1 + h) = (R/R)(a/a), 
where R (t) is the radius of the spherical region; similarly we may approximate £ by the 
density contrast 5. Both these approximations are valid for £ ^ 1 and should be interpreted 
as some average value. From the STH model it is easy to show that (with £ = <5), 

1 | <F = (/i-sin/i) 2 8(0- l) 3 

(1-cos^) 3 Q3 [cos -i ( 2_ l) _2 v /^- 



and 



l + fr= ( ^ Sin//)si 7 , ^f"^, ,, (17) 
(1-cos^) 2 ft [cos-* (£ - 1) - gVST^T] 



where // is the phase angle describing the evolution of the STH. It follows that 

(l + /i) 2 fl sin 2 n 



(18) 



1+f ~~ 2 (1 - cos^)' 

This is valid for all values of equations (16) and (17) can be analytically continued for 
O < 1. A nonlinear structure which is undergoing free-fall collapse, we will have \x w 2n — e 
with small e and £ ^> 1, h 2 ^> 1. In this limit equation (18)becomes 

h = Kf 2 ft 1/2 - (19) 

The subscript "eff" is a reminder that we should use an effective O which exists around 
most collapsing structures at the epoch of interest. [We are assuming that the shperical 
top hat region can be modelled as being located inside a larger region, which - when 
smoothed over the relevant scale - behaves like a universe with an effective ft e ff-] Though 
it is difficult to model this quantity, two features are apparent: (1) At a given epoch, 
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h oc in the nonlinear phase. (2) As clustering produces underdense regions more and 
more collapsing structures will find themselves in Q e f / < 1 regions. As time goes on Q e f f 
should decrease, and h should decrease at a given £. 

Figure 3d shows that both these expectations are borne out. The dashed lines at the 
nonlinear end has the slope of 0.5 which fits quite well to the simulation results. At a given 
£, the h decreases as time goes on; at z « 0, we find Q e ff ~ 0.3. 

In the case of Figures (3a,b) the key point to note is the lack of stable clustering. 
In both the cases, h does not approach unity asymptotically. This fact is more easily 
discernible using h than from mean correlation function. The value of h at a given £ 
decreases as we add more small scale power, as expected. Note that lack of stable clustering 
invalidates the argument based on (12). Also note that these results disagree with the 
results of simulations discussed in Mo et al. (1995). They claim £(a,;r) oc [£i(a, Z)] 3 / 2 
asymptotically, which requires h = 1 for large £. Our simulations show that ft ^ 1 
asymptotically. 

The models described above all have = 1. To check the dependence on the 
background scale, we considered the behavior of h[£\ in the following cases: (i) a flat 
CDM + A model with Vic dm = 0.4 and Oa = 0.6 (model e) and (ii) an open model with 
O = 0.1 and = Ak~ l (model f ). These models have specific background scales and 
we find that it is reflected in h[£] [see Figures (3e,f)]. The deviation in the linear end is 
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understandable since 

s -^=§(s) « 2o) 

where b (t) is the growing solution to the linear perturbation equation. For reasonable fi, 
we can approximate (ab/ab) by [fi (z)] ' 6 . This formula is used usually at z = 0; it is a 
reasonable approximation at moderate z, if the instantaneous value of O is used. So, in 
the linear end h oc [Cl (z)] ' 6 at a given £ and this is indeed seen in figure 3 (e, f). At the 
nonlinear end, h seems to approach unity within the limits of our resolution. This suggests 
that models with O < 1 obey the hypothesis of stable clustering to a better approximation 
than models with = 1. This is understandable since structures freeze out in the < 1 
model around z ~ O -1 — 2; in contrast, significant merging takes place in O = 1 universe 
at all epochs. In the intermediate scales, no clear pattern emerges in these models. 

4-3 Relation between £ and £l 

The original ansatz can also be presented as a relationship between £(a,x) and 
£z,(a, Z). In the case of strict universality, these two approaches are related by equation 
(9). However, when deviations from universality are present, it is worthwhile to examine 
both independently. Figures (4a,b,c,d,e,f) show the results for different cases. These 
figures relate the nonlinear correlation function at a scale x to linear correlation function 
at a scale / where / = x(l + 1 / 3 . The relation between / and x is independent of the 
value of O and can be derived by the method of characteristics used by NP. The broken 
lines in these figures have slopes 1 and 1.5 as suggested by equation (10). We find that 
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the deviations from our ansatz, seen in h, give rise to corresponding deviations in these 
figures. A spectral dependence is evident in Figures 4a and 4b, as well as in the slight shift 
for the CDM case, Figure 4c .This dependence seems to be in the same sense as reported 
by Mo et al, 1995. The deviations from stable clustering follows the pattern discussed 
earlier, except for Figure 4f. When O ^ 1, we have £l = b 2 (t)x~^ n+3 ^ and £ = a 3 F(ax); 
in that case, stable clustering assumption does not translate into a simple scaling of the 
form £ oc £^/ 2 . By and large, these figures reinforce the conclusions arrived at earlier. 

4-4 Overall Conclusions 

1) As expected, (but contrary to some earlier work), we find that if the input Pk 
has a power law dependence on k, then the only scale that is relevant is the scale at 
which mass fluctuations are nonlinear, x^l, which induces a break in the input (linear) 
power spectrum. Under these circumstances all quantities that we have measured show a 
self-similar evolution, i.e., they are only functions of scaled radius x/xnl- 

2) The "stable clustering" assumption is not true in general, presumably because, in 
the case of O = 1, there is steady merging of structures so that isolated bound objects do 
not survive intact, violating the key ansatz behind the stable clustering hypothesis. When 
O < 1, since structure freezes out at low redshift, the stable clustering model is a better 
approximation than if O = 1. 

3) The relation between the dimensionless pair velocity h(a,x) = —v pa i r /Hr and 
the mean correlation function £(a,x) is not universal because of the failure of the stable 
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clustering assumption although within the range of cosmologically plausible models the 
relation h(£) shows sufficiently uniform behavior to be a useful quantity. 
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FIGURE CAPTIONS 



Fig. 1- Figures (la,b,c,d,e,f) show the two-point correlation function as function 
of r/r n i (where r n i is the scale corresponding to £ = 1) at three 
different redshifts [z = 0, 0.5, 1 with (solid, dotted, dashed) curves, 
respectively; this order will be maintained in subsequent figures] for 
Models (a,b,c,d,e,f) listed in Table 1. Also shown as two dashed lines 
in each figure are asymptotic slopes at the extreme nonlinear and linear 
ends. Note that figure 1 (a,b), which are for powerlaw spectra in O = 1 
universe, show the expected self similar behaviour while the other figures 
do not. 

Fig. 2- Figures (2a,b,c,d,e,f) show the ratio f3d(«, x)/[ax non u nea r («)] between 
the 3-dimensional pairwise velocity dispersion to our velocity unit, as a 
function of r/r n i at three different redshifts (z = 0,0.5,1) for Models 
(a,b,c,d,e,f). Figure 2 (a,b), which are for powerlaw spectra in O = 1 
universe, show the expected self similar behaviour while the other figures 
do not. 

Fig. 3- Figures (3a,b,c,d,e,f) shows h as a function of £ at three different 
redshifts [z = 0,0.5,1 with (solid dots, open circles, open squares), 
respectively; this order will be maintained in subsequent figures] 
for Models (a,b,c,d,e,f). The errorbars are shown for la statistical 
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uncertainties. In the nonlinear end the h is lower for a given £ if there is 
more small scale power. The value of h in the asymptotic limit depends 
on the spectrum and does not always go to unity; this shows that stable 
clustering is not a valid assumption in general. 
Fig. 4- Figures (4a,b,c,d,e,f) show £ as a function of £l at three different redshifts 
(z = 0,0.5, 1) for Models (a,b,c,d,e,f). 
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Table 1. Summary of the computed models 



Run 


Model 


Code 





A 


L box h~ 1 Mpc 


Iresh 1 MpC 


&8 


a 


Pk = Ak~ 2 


PM 


1 





128 


0.18 


1.05 


b 


Pk = Ak- 1 


P 3 M 


1 





128 


0.025 


1.05 


c 


SCDM 


P 3 M 


1 





128 


0.025 


1.05 


d 


HDM 


PM 


1 





128 


0.18 


1.05 


e 


LCDM 


P 3 M 


0.4 


0.6 


128 


0.025 


0.79 


f 


Open P k = Ak~ l 


PM 


0.1 





128 


0.18 


1.05 
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